Applied Probability Trust (23 April 2012) 



A WAVELET-BASED APPROXIMATION OF FRACTIONAL BROW- 
NIAN MOTION WITH A PARALLEL ALGORITHM 

^ DAWEI HONG,* Rutgers University 

^ SHUSHUANG MAN," Southwest Minnesota State University 

JEAN-CAMILLE BIRGET,* Rutgers University 
DESMOND S. LUN,* Rutgers University 

Abstract 

We construct a wavelet-based expansion to approximate fractional Brownian 
motion of Hurst index H £ (0, 1). For practical implementations, the expansion 
converges almost surely and uniformly in discrete time t € [0, 1]. We prove that 
the convergence rate is optimal. We also show that the approximation can be 
implemented by a fast parallel algorithm. 

Keywords: Fractional Brownian motion; wavelet expansion of stochastic inte- 
gral; uniform approximation; convergence rate 
2010 Mathematics Subject Classification: Primary 60G22 

Secondary 65T60; 65Y05 



* Postal address: Center for Computational and Integrative Biology, Department of Computer Science, 

Rutgers University, Camden, NJ 08102, USA; Email address: dhong@camden.rutgers.edu 

** Postal address: Department of Mathematics and Computer Science, Southwest Minnesota State 

University, Marshall, MN 56258, USA 



1 



2 



D. Hong, S. Man, J-C. Birget, D. S. Lun 



1. Introduction 

A fractional Brownian motion (fBm) (-Bf^'')tG[o,T] of Hurst index H € (0,1) is a 
centered Gaussian process with covariance Elslf' B^^^] = (l/2)(t^^ + — \ti ~ 
<2p^) for all t\^ti G [0,r]. A standard Brownian motion (Bm) {Bt)te[o.T] is the 
special case of iJ = 1/2. There are enormous applications of fBm in engineering and 
sciences; see |3] and references therein. The study on approximation of fBm has been 
active since the 1970s. A major focus is to find approximations of fBm that converge 
in law; for example, see [5J (THl [H [H [H] and references therein. However, practical 
implementations often require almost sure and uniform approximations of fBm. Such 
an approximation works as follows. Let {Bj. )te[o,i] be a fBm of some H e (0,1). 
Then, with respect to the probability space where {Bj. )te[o,i] is defined, the following 
event occurs with probability 1. For a sample path of {Bj )te[o,i] there is a sequence 
of functions of i G [0, 1] produced by the approximation which uniformly converges 
to the sample path; conversely, a sequence of functions of t G [0, 1] produced by the 
approximation uniformly converges to a sample path of {Bj. )tg[o,i]- Meyer, Sellan 
and Taqqu |17) obtained a wavelet series expansion of fBm which provides an almost 
sure and uniform approximation of fBm of iJ e (0, 1). The result also brings deep 
insights into spectral properties of fBm. The authors showed the existence of optimal 
wavelet expansion of fBm. However no convergence rates were given. Kiihn and Linde 
|12) proved that the optimal convergence rate that any series expansion of fBm may 
reach is 0{N~" y^log N) if the expansion converges uniformly over time interval [0, 1]. 
Dzhaparidze and van Zanten [7] constructed an explicit series expansion of fBm of H 
€ (0, 1). With probability 1 this series expansion converges absolutely and uniformly 
over time interval [0, 1]; and soon after the authors proved that the convergence rate 
is optimal [8]. 

The above results will have a long lasting impact on the study of fBm. Nonetheless, 
there are two questions. The series expansion of fBm constructed in 7 is based on 
{sm{xnt) /xn)n>i SLTid (cos(?;„t)/j/„)„>i where x„ and y„ are the positive zeros of the 
Bessel functions J-h and Ji^h respectively. The wavelet series expansion of fBm in 
[17) uses the Haar wavelet; and in the remark on Theorem 2, the authors pointed out a 
technical difficulty to use the Mandelbrot - van Ness stochastic integral representation 
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of fBm 16 for wavelet series expansion of fBm. In practical implementations, the 
Haar wavelet is easier to compute than positive zeros of Bessel functions. Moreover, 
the simple form of the Mandelbrot - van Ness representation is likely to yield a fast 
algorithm. Thus, a question is if, using the Mandelbrot - van Ness representation, we 
can find a Haar wavelet series expansion of fBm of H e (0, 1), and prove that the 
expansion provides an almost sure and uniform approximation with the above optimal 
convergence rate. 

In the present article, we construct a wavelet series expansion of fBm of iJ G 
(0, 1) that meets the requirements mentioned above. Our method is to apply Levy's 
equivalence theorem (see Theorem 2.6.2 in [13J to a Haar wavelet expansion of fBm 
obtained based on the Mandelbrot - van Ness representation, and then to carefully 
evaluate the wavelet coefficients. We focus on discrete time t E Q O [0, 1] where Q 
stands for the set of rational numbers. On one hand, the set of rational numbers covers 
all needs by practical implementations. On the other hand, by the Holder continuity 
of fBm [T] with probability 1 

lim sup — = Ch for a constant Ch > depending only on H 

and hence, our result can be extended from the case oi t G Q n [0, 1] to the case of i G 
[0, 1]. Moreover, we show that our constructed wavelet series expansion of fBm yields 
a fast parallel algorithm for almost sure and uniform approximation of fBm of _ff G 
(0,1). 

Another question is if the constant behind the big O in the optimal convergence rate 
is related to Hurst index H . When a convergence rate for a series expansion of fBm 
is written as 0{N^^ y/\og N), the fact that fBm cannot be defined for i7 = 1 cannot 
be seen. For our wavelet expansion of fBm, we show that the constant behind the big 
O is of the form Cy/q{H{l — H))^^/^ where C is absolute and q > 1 is a parameter 
in the large deviation bound for the almost sure convergence; see Theorem 15.11 - the 
main result in this article. 

Note that almost sure and uniform approximations of fBm were investigated by 
other approaches. One approach is by moving averages of simple random walks with 
convergence rate 0{N-^°s4 mi„{ff-i/4, 1/4} \ogN) only for H G (1/4, 1) [H]. Another 
approach is by transport processes with convergence rate 0(iV^(^/^~'^)(logiV)^/^), \H— 
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1/2| < P < 1/2, for H E (0,1) in 110 . The fastest convergence rate for almost sure 
and uniform approximations of fBm oi H G (0,1) appears to be the optimal rate by 
series expansions of fBm. 

The rest of the present article is organized as follows. Section [5] is for preliminaries. 
In sections [31 2] and [5] we construct and prove a uniform approximation (in discrete 
time) of fBm of Hurst index H e (0, 1). In section [HI we present a parallel algorithm 
for the constructed approximation. 



2. Preliminaries 

The Mandelbrot - van Ness stochastic integral representation of fBm [TB] is 

Bi"^ = Ch f i{t - s)l-'" {-s)^-''') dB,, t e [0, 1] 

Here H is the Hurst index taking values in (0,1) and Ch = (r(-ff + 1/2))^^, the 
reciprocal of the Gamma function ai H + 1/2; sec 1.2 in 3 for results from further 
studies on the above representation of fBm. In what follows, we denote the underlying 
probability space for the above representation of fBm by (17,7^, Pr) where is a 
standard Brownian filtration. 

Our construction of a uniform approximation of fBm is based on a rewriting of the 
Mandelbrot - van Ness stochastic integral representation 

Bf^ = h{t,H) + h{t,H)+h{t,H), te[0,l] (1) 

where 

h{t,H) = CH f\t-s)"~''HB, 
Jo 

hit, H) = Ch 1° ((i - s)^-^'^ - {-s)"-'/^) dB, 

h{t,H) = Ch r Ut - sf-'/^ - i-sf-'/') dB, 

Let {(t>n)n>o be a complete orthonormal basis for L^[a, 6]. For / e i^[a, 6], / — 
(/' '^n) 't'n in Li^ [a, • Take Wiener integration on both sides of the above equality, 
and informally interchange the order of integration and summation. We have 

„b oo „b 

/ /(s)dB, = ^ (/,</>„) / <tyn{s)dBs (2) 
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We can check that the right side of ^ converges to the left side in E[| • p]. However, by 
Levy's equivalence theorem we ftirther have the following (see Theorem 2.6.2 in [T5]): 

Theorem 2.1. The right side of ^ converges almost surely, i.e., equality ^ holds 
with probability 1 . □ 

The Haar wavelet on [0, 1] is defined as follows. Let 



nis) = < 



1 if0<s<l/2 
-1 ifl/2<s<l 
otherwise. 



For n = 23 + k with j > and < < 2^ define Hnis) = 2i/'^'H{2is - k) and 
Hois) = 1. The sequence {'Hn)n>a is the Haar wavelet on [0,1], which constitutes 
a complete orthonormal basis for L^[0, 1]. In a similar way, we can define the Haar 
wavelet on any given interval [a, fe] C M to constitute a complete orthonormal basis for 
L^[a,b] (see [4]). 

3. Approximation of /i (t, jy) 

In this and the next section, we construct and prove approximations of Ii{t, H) and 
/2(i, -ff), respectively. The two approximations are uniform over time t £ [0, 1] H Q. In 
the meantime, we investigate how the convergence rates of these approximations are 
degrading when Hurst index iJ — > 1_ or 0+. 

Consider a family of functions in L^[0, 1], /j^"* with a parameter t e (0, 1] n Q. For 
every t g (0, 1] n Q 



{t~s)"-^/^ iise[0,t) 

otherwise. 



By Theorem 12. II we have 



Pr|(^'/**'^(5)rfB.) (^) = ^f2yi'\nn) j\r.{s)dB. 



(^) = 1 (3) 
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for each t G (0, 1] fl Q and as a consequence 

Pr fl I f rfi'\s)dB,) M = (f^if^'KUn) fun{s)dBA = 1. 

tG(0,llnQ I ^-^O ^ \n=0 -^0 / J 



E(0,l]nQ 
We define for all TV > 1 

W^{t,H,N) = I 



Ch T.n=M^\T^n)C^n^ for t € (0, 1] H Q 

(4) 

for f = 



Here 6^^ = Hn{s)dB„ n = 0,l,...,N, are i.i.d. Gaussian random variables with 

mean and variance 1. 

In what follows we use two conventions: n e Z+ is said to be at level j if n = 2^ + 

k with j > and < fc < 2-'; and the interval is meant to be when 

fe+i _ 1 

Lemma 3.1. There is an absolute constant Di > such that for every t G (0, 1] fl Q 
and for all N > 1 



Proof. For t G (0, 1] fl Q, at each level j = 0, 1, . . ., we partition the set 

{n = 2^ + A;:fc = 0,l,...,2^-l} 



into three subsets: Qi{j,t) consisting of all n (if any) such that t > Q2{j,t) 
consisting of one n such that t G -^ti); and Gaij, t) consisting of all n (if any) such 
that t < ^. 

Consider a fixed j. By the definition of /j^' we have for n e Qs{j,t) 



fi'\Hn)=0 (5) 

Let us denote by 2^ + kf^ the unique n G G2{j,t) such that t G [%^, -^2r~)- Then 
we have n = 2^ + ktj and 



A wavelet-based approximation of fBm 



7 



By the above equality we see that 



H-l/2 



/ 2kt 



H-1/2 



2JH 



ds 



By this inequahty, with calculation we have for n Cz G2 {j, t) 



fr,nn) < 2-2J^^ 2-(2ff+i) (if + 1/2) 



For each n £ Qi{j,t), we have n = 2^ + k with k < ktj and 



/i^\H„)=2^-/2 



2fc + 2 



j 2k I 2k + l 



2^/2 



H + 1/2 



t - 



2fc 



2JH 



H+1/2 



- u- 



2fc + 1 



2fc + l\«+^/^ 



t 



2J+1 

2fc + 2\^+^/n' 



2J+1 J \ 2J+1 

To facilitate our argument, we introduce a function w oi h 



(6) 



(7) 



= 5(2:0 + h)+ g{xo -h)~ 2g{xo) 



where 



ff(') ^ (•)^^"'^''^ and xo — t 



2k + 1 
2J+1 



Then, letting ft, = 27Tr' we rewrite ([7]) as 



Consider the following Taylor's expansion 



2J/2 



H + 1/2 



w{h) 



w{h) = w(0) + 



w'(0) w"(6>fe) 
1! 2! 



ft^ (for some < 6* < 1) 



w"{eh) 
2! 



ft,2 (since w(0) = ^'(0) = 0). 



Hence 



w{h) = 



'{Oh) 



2! 



= 2-2(j+i) 



(i7 + l/2)(i?-l/2). 



2J^ 



2J+1 



(8) 
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Now consider those n e Gi{j,t) with n = 2^ + k and k + 2 < ktj. By the above 
equahty and ([8]) we have 



< 2^/22-2(j+i)|j^_ i/2| [t- 



2k + l + e\"^^^^ 



since < 6* < 1 and < H < 1. Furthermore, since t e ^i^), by the above 



equahty we have 



It : rLn 



H-3/2 



_ |g-l/2| ^_,g 



(fc:}-(fc + l)) 



H-Z/2 



Thus, for 11 G Gi{j,t) with n = 2-' + and fc + 2 < /ctj-, we have 

2H-3 l/2|2 



16 



(9) 



There is one and only one (^f^^\T-Lj^ with n G Gi{j,t) which is not included in Q, 
namely the case of n = 2^ + kt.j — 1. However, in this case we have 



/«,H„U2^-/2 



and hence 



2J 



.2-(2ff+ib- ^2-2j"^(i7 + l/2)- 



(10) 



(-ff + 1/2)2 

Now, putting (O, ([S]), (jni and together, we have that there is an absolute 
constant £>! > such that at any level j 



E 

{n at level j} 



Without loss of generality, we can rewrite the above inequality as 



E 

{n at level j} 



/t I 



- 1 -i? 
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Therefore we have 



E (/.^"^^"> ^ E E 

n=Af+l j=Llog2 N\ {11 at level j"} 



f(l) a/ 



< 



^ 1-r 



l-H 



1 



H 



1 -i? 



1 - 2 



-2ff 



Since HmH-i.o+(l ^ 2^^^)/H = 2 log 2, there is an absolute constant G > such that 
1/(1 — 2^2^) < G/H for all iJ e (0, 1). Hence we can rewrite the above inequality as 
claimed by the lemma. □ 



Lemma 3.2. For any given H e (0, 1) and q > 1, we have for all N > 1 



pj »p |/.(,H)-HM<.H.iV)|> ""f^'-^ 



< 



1 



where Di is the absolute constant used in Lemma lS.l 



Proof. By definition Ii{0,H) = = Wi{0,H,N). So, we focus on the case of t e 
(0, 1] n Q. By ([U, ([3]) and its consequence we have 



Pr fl {{h{t,H)-Wiit,H,N))iuj) 
te(o,i]nQ 



00 „i 



(11) 



n=N+l 



Here X]^Ar+i(/t ^n) Jq 'Hn{s)dBs is a Gaussian random variable with mean and 
variance X]^A^+l(/^^^^n)^■ We denote this variance by (j\{t,H,N). 
For any given H e (0, 1) and q > 1, we have 



{00 „i 
E / (ft'\nn)nnis)dB, 



> 



y/2Diq\ogN 



7]- / ^2Jii;1oeJV 
W^^/f (l-H) 

^ 2a,{t,H,N) ^^^ 



exp 



2a2(t,i/,7V) 



N"y/H{l - H) 

du 



/TT 



u 



VSi'igiogTV \^/2ai(t,H,N) 

nH ^ H{l-H) 



exp 



<Ji{t,H,N) ^ N"^H{l-H) 



X exp 



V2ai{t,H,N) 
DiqlogN 



du 



al{t,H, N)N^"H{1-H) 



< 



1 



V27rglogiV 7V« 



by Lemma 13.1 
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By inequality and since 1/ \/2ti log < 1 for all > 1, we have completed the 
proof of the Lemma. □ 

A remark for the above Lemma is as follows. By Lemma and the Borel-Cantelli 
Lemma, we can see that when N ^ oo, Wi{t, H, N) converges to Ii{t, H) almost surely 
and uniformly in < G [0, 1] n Q. In the meantime we also see that due to the factor 



\/ yjH{\ — H) the convergence rate may degrade as the Hurst index H ^ 1_ or 0+. 
The result in [12^ shows that the optimal convergence rate is 0{N^^ y/log N). Hence 
the above factor causes a degradation of the convergence rate when H ^ 1_ or 0+. 



4. Approximation of /2(t: -f^) 

The construction of an approximation of I2{t^ H) follows a procedure similar to that 
in the previous section. Consider the Haar wavelet ('H„)„>o on [—1,0]. Define a family 
of functions in 1,0], /j^"* with a parameter t e [0, 1] n Q: 



By Theorem 12. II we have 



(i_s)ff-l/2_(_s)H-l/2 ifse[_l,0), 

otherwise. 



Pr I /f (io) = (^fy.fi'\nn) £nnis)dB}j = 1 (12) 



for each i e [0, 1] n Q, and as a consequence 



Pr fl Iff" fi^\s)dB^ (oj) = lf^{fP,Hn) f Hn{s)dBs 
te[o,i]nQ I ^ ^ \n=o -^^1 



(lo) = 1. 



We define for ah A^ > 1, 



, Ch En-o if , nn)clt> for t e [0, 1] n Q 

W2{t,H,N) = { * (13) 

for i = 



Here ci? = /^^H„(s)dB„ n = 0, 1, . . ., A^, are i.i.d. Gaussian random variables 
ar 

sequence (£« ^)n>o used in the definition of Wi{t, H, N). 



(2) 

with mean and variance 1. Notice that the sequence {Cn )n>o is independent of the 



A wavelet-based approximation of fBm 



11 



Lemma 4.1. There is an absolute constant D2 > such that for every t G [0, 1] H 1 
and for all N > 1 

2 Do 



Proof. For each t e [0, 1] n Q 

n=N+l 



/ 00 

<2 E E ((--)^-^/^^« 



(14) 



\ii=JV+l n=N+\ / 

Here ((i — s)^^^^"^ and {{—s)^^'^/'^,'Hn) are Riemann integrals. Thus, instead 
of ((i — s)-^^^/^ jHn) and {{—s)-^^^^^,'Hn), by changing variables we can estimate the 
right side of the above inequality, by using {{t + s)-^^^/^ jHn) and (s^^^/^ Below, 
we estimate Er=Ar+i((i + s)^"^^^: ^«)^ and Er=A'+i(s^"^^^''^n)^ separately. 

We estimate the first term. For t e (0, 1] n Q, at each level j, we have for each n = 
2^ + fc, fc = 0, . . ., 2^ - 1, 



2j72 



L 2J- 
2^/2 

i/ + l/2 



+ sf-'/-'ds - r-' {t+ s)"-'/^ds 



t + 



2k + 1 



2J+1 



t + 



2J+1 



t + 



2fc + l 

- \t + 

2fc + 1 



2J^ 



2fc 

H+l/2' 



H+l/2\ 



(15) 



To facilitate our argument we introduce a new function w of h which is a revised version 
of the function w in the proof of Lemma 13.11 



where 



w{h) = 2g{xo) - g{xo + h) - g{xo - h) 



9i-) = i-) 



_ , -.H+1/2 



and xq = t 



2k + 1 
2J+1 ■ 



Then, letting h — 2jTr, we rewrite ([T5l) as 



((^ + ,s)^-l/^H„) 



2^/2 



w(/i) 
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Then, similarly to the proof of Lemma |3.1[ we consider the following Taylor expansion: 



w{h) = w(0) + 



w'(0) w"{eh) 



2! 



1! 2! 

(since w{0) = w'{0) = 0) 



h'^ (for some < 6* < 1) 



Hence we have 



w{h) 



t + 



^ 2 w''{0h) ^ 2f,+i) (-g + V2)(g-l/2) 
2! 2 



2J+1 



+ \t + 



2J+1 



and consequently 



((t + < V/^2-W)\H - 1/2| (^t + ^^i^) 



(16) 



since < 61 < 1 and < iJ < 1. Now, as in the proof of Lemma we denote by 
S-* + fctj- the unique n such that t e ["Ir", Then, there are two and only two 

cases. 

Case 1: ktj > 1. By (fT6|) we have 



(^ + .)^-l/^H,: 



H-3/2 



2^+1 2J+1 / 

< 2-'^" \H ~ 1/2| 2-(^+i/2) {k + 1)^-3/2 
Thus, there is an absolute constant I?2,i > such that 

2 °° / 1 

((^ + .)^-^/^H„) <^2.i2-^^-^E(7 

{n at level j} 



e=i 



3-2H 



(17) 



Case 2: kt.j = 0. Using (fT5|) we have 

((t + .)«-V2^^^^^ 



2^/2 



{t + s)"-^'^ds~ r' it + sf-'/^ds 



2^/2 



i? + 1/2 
2 



X H+1/2 



2J+1 

H+l/2 ^ -j^ X _ff+l/2 



/ 2 / 1 y 
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< .Jl^ f 1] (18) 



and hence 



For n = 2^ + k with = 1, . . ., 2^ - 1, by we have 

(19) 

< ^2-20+1) _ since - 1/4| < 1 for e (0, 1). 

Putting (fT8|) and (|19p together, we have 

Y: {{t + sf-'/^Ur) <D,,,2-^^"Y.[l 

{n at level j} f=l ^ 

Here, without loss of generality we can let -D2,i be the same absolute constant as in 

m- 

With an argument similar to the above for ((t + s)-^^^^^ jHn), we have that there 
is an absolute constant D2 2 > such that 



(20) 



2 00 ^ ^ 3-2H 

{n at level j} 1=1 ^ ^ 

The lemma follows from putting ([14]), ([20]) and ([21]) together. □ 
Lemma 4.2. For an?/ given H £ (0, 1) and g > 1, we have all N > 1 



(21) 



where D2 is the absolute constant used in Lemma \4.1\ 
Proof. By (fT3|). (fT2|) and its consequence we have 

Pr fl {ihit,H)-W2{t,H,N))iLo) 

te[o,i]nQ 

(-0 



Cff ^ (/f\H„) / H„(s)di?,(c^)} = l 

n=Af+l -^^1 



Here, J2'^=N+iift^^ ^'^n) j'^i'Hn{s)dBs is a Gaussian random variable with mean 
and variance X]^JV+l(/^^^^")^■ -By Lemma 14711 following the same procedure as in 
the proof of Lemma 13.21 we complete a proof of the lemma. □ 

The result of Lemma 14.21 has an expression similar to that of Lemma 13.21 Thus the 
remark made there holds here. 
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5. Approximation of fBm 

For H e (0, 1), given < ti < . . . < U < . . . < U < 1 with all U e Q, find a 
sequence {W{t, H, A^))jv>i of functions defined on i e [0, 1] such that with probability 

( M ) 

1 the above sequence converges to a sample path of (Bl )fg[o,i] uniformly over t — ti, 
i = 1, . . ., i. For a given N, W{t, H, N)) is called the iVth step of a uniform in discrete 
time approximation of {Bl )te[o.i]: simply the iVth approximation of {Bf )tg[o,i] 
when all U are regarded as predetermined. With this in mind, let us take a close look 
at {Bf^^)t^[o,i] represented by H]). 

h(t,H) is the process (Ch C(t - s)"-^/'^dB^ . At each time t e [0,1] the 

V " ^te[o,i] 

process is defined by an Ito integral Ch /q — s)^^^/^dBs which itself is a process on 
[0,i]. We considered 

Ch f {t^s)"-''^dB^ (22) 
"'o /te[o,i]nQ 

and defined Wi{t,H,N) that converges to almost surely and uniformly in t e 
Qn [0,1]. Here Wi{t,H,N) is defined via c';l'^ = J^nn{s)dBs, n = 0, 1, . . ., and 
hence, it is defined on (i7,J^, Pr). 

l2{t,H) is the process f Ch f° ((t - s)"-^/^ - (s)"-^/^) dB^) . This pro- 

V ^ 'I tG[0,l] 

cess is indexed by time in the usual way, i.e., at each time t e [0, 1] the process is defined 
by the random variable Ch jli {{t — s)^~^^^ — (— s)^^^/^) dBg (not an Ito integral). 
Thus 

Ch f ({t s)"-'/^ i-sr-"^) dB\ (23) 
"'-1 / ie[o,i]nQ 

is a process defined on (ri, J^, Pr). Because the integrand tends to infinity as s — > 0_, 
in the case oi H G (0, 1/2) we defined W2{t, H, N) that converges to ([231) almost surely 
and uniformly in i G Qn [0, 1]. Here W2(t, H, N) was defined via c'h'' = /"^ 'H„(s)dB^, 
n = 0, 1, . . ., and hence, it is defined on {ft,T, Pr). 

Isit, H) is the process f Ch {{t - s)"-^'^ - {-s)"-^''^) dB,) . This pro- 

V °° / 46 [0,1] 

cess is indexed by time in the usual way, i.e., at each time t G [0,1] the process is 
defined by the random variable Ch JI^ {{t — s)^^^/^ — (— s)-^^^/^) dBg (not an Ito 
integral). Thus 

(ch r (it s)"-'/^ {-sr-'/A dB^l (24) 

V J-oo ^ ' J tg[o,l]nQ 
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is a process defined on {fl,J-, Pr). Since the integrand is well defined in s G (—00, — 1], 
we have a direct simulation (not an approximation) of (|24p . Combining this sim- 
ulation of (p4|) with Wi{t, H, N) and W2{t, H, N), we obtain an approximation of 
(-Bj )tG[o.i]nQ- Below we carry out this idea. 

For a fixed t G [0, l]nQ, Ch JI^ {{t - s)"-^/^ - {-s)"-^/^) dB, generates a Gaus- 
sian random variable X{t) with mean and variance 

J —00 

as follows. For almost every cu ^ 

X{t){u) = Ch I \{t- s)"-'/' - hs)"-'/A dB,{u:) (25) 

J —00 

By definition X{t) is independent of Cn^ and Ln \ n = 0, 1, . . ., in the definition of 
Wi{t,H,N) and W2{t,H,N). By <^ we have 

Pr n {ch r\{t~sr-'/^-{-sr-'l^)dB,{u.) = 
te[o,i]nQ ^ °° 

Ch (/"J ((t- .)^-V2 _ (-,s)^-V.)^rf,y^'^ = 1 

where £ is a standard Gaussian random variable independent of c\}^ and L^n \ n — 0, 

1, Hence, given ti e [0, 1] n Q, i = 1, . . ., i^, we can directly simulate l3{ti,H) as 

follows. Generate a sample C{uj) of a standard Gaussian random variable independent 
of Cn^ and n = 0, 1, . . .; then for i = 1, . . ., i, let 



h{U,H){uj) = 




Now, in (51, Pr) we define 

W(t , if, iV) = Wi (t , ii, iV) VF2 (t , H, N) +h{t, H) for i e [0, 1] n Q 



Theorem 5.1. There is an absolute constant C > such that, for any given H G 
(0, 1) and q > 1, we have for all N > 1 



Pr. 



sup 

te[o,i]n 



B. 



W{t,H,N) 



> 



y/H{l - H) 



< 
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Proof. Using 2N + 3 independent samples from the standard Gaussian distribution, 
we can approximate a path in a fBm {Bl )te[a.i\ foi' given ti G [0, 1] H (Q, i = 1, 
as fohows. We use the first + f samples to obtain I instances of Wi{t, H, N) 
ai t = ti according to and use the next + 1 samples to obtain t instances of 
W2{t, H, N) a.t t — ti according to p^ . By definition, these 2£ instances approximate a 

( H ) 

common path in (B^ )tiE[o.i]- We use the last sample to obtain i instances of Isit, H) 
aX t — ti according to ((26|) . By definition, the obtained t instances belong to the path 
in {B\ )fg[o,i] which is approximated by using the first 2N + 2 samples as described 
above. Then the theorem follows from Lemma [3.21 and [4.21 and the fact that l 

can be any number in N. 

Before leaving the proof, we note that the variance in (PHl) is bounded for t e 
[0, 1] n Q. Indeed, 



Jo 

<\H- l/2|s^-3/2 f du = t\H - l/2|s^-3/2 
Jo 

and hence, 

d p ({t~ s)^-i/2 _ i-s)"-^'^) ' ds = Cl (it + s)"-^'^ - ' ds 

„2H-3 j„ _ Cjjt^{H - 1/2)2 



/OO 
s'^-^ds 



□ 



2{l-H) 

6. A parallel algorithm for the approximation 

In order to specify a concrete approximation of fBm, an input of the following type 
is given: 

• A sequence of ^ time instances Q < ti < . . . < ti < . . . < t^ < 1\ 

• an approximation bound a > 0; this is a probable upper bound on the distance 
between fBm and the approximation; 

• a probability bound p > such that 1 — p is a lower bound on the probability of 
correctness of the given approximation bound a. 

• Optionally, the Hurst index H may be taken as an input, or may be viewed as fixed. 

The algorithm has two parts. First, a preprocessing step uses a and p (as well as 
H) to find the corresponding parameters N and g, in order to obtain the deviation 
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bound of Theorem 15.11 Second, a probabilistic parallel algorithm takes as inputs A^, 
q, H, and sample times ti (i ~ 1, ...,£); it outputs a sequence of numbers W{ti, H, N) 
{i = l,...,l) such that Pr{supi<,<^ \B[f^ - W{U,H,N)\ > a] < p. 

Part 1: Preprocessing step 

On input a, p (and H) we find N and q by solving the system of two equations, 
obtained from Theorem 15.11 



(A) a = ^ 

(P) p = 7m 

This system is always solvable for p e ]0, 1] and small enough a > 0, which can be 
seen as follows. We eliminate N and derive an equation with unknown q: 

_ C ^(l/g)-log(3/(p7g)) 



The numerator y^{l/q) ■ log(3/(pY^)) decreases strictly from •\/log(3/p) (when 
g = 1) to (as g — 7^ +oo). The denominator {i/ {p^/q))^^'^ is bounded for q e [1, +oo] 
(it tends to 1 as g — > +oo), and is also bounded away from for any given value oi p 
and H . This implies that equation (Q) is solvable in q for any p G ]0, 1], G ]0, 1[ and 
a g ]0, aijiax[ (for some Omax = CLiaa.x{H,p) > 0). Once a value for q has been found, 
equation (P) yields a value for N (which we round up to the nearest integer). 

For the complexity analysis an upper bound on in terms of a and p is needed. 

Proposition 6.1. Let N be the solution of equations (A) and (P). Then \ogN < 
■ (llog(pa)l + log -7^^). Thus for fixed H, logN < 0(| log(pa)|). 
Moreover, logN > -g^-(|loga| +log— ^===y), hence for fixed H , \ogN > f2(|loga|). 



Proof. We have > -y/log N, for all £ > 0, iV > e. We will also need H > e, so we 



pick e = H/2. Equation (A) then implies < (- 



Replacing y/g — (from equation (P)), we obtain 

/ 3c 1 .2/H 1 . 3c 1 1 

N < { , — • — ) ttjt; < ( . — — I — £7777 smce 1 < g 



Hence, 



J^l + H/2 ^ ( 3c . J_^2/-f^ 

^ y/Hil-H) pa 



The first inequality then follows immediately. To obtain the last inequality we use 
< ^/TogN and 1 < g in equation (A). T 
the proposition follows immediately. □ 



1 < \/\ogN and 1 < g in equation (A). This yields > ^^^'^ ■ i, from which 
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Part 2: Parallel algorithm to find approximants at sample times 

In the following parallel algorithm, only step 1 is probabilistic. 

In this part, only N and £ are used as complexity parameters. Since an approx- 
imation bound a is given, the numerical representations and calculations require an 
absolute precision < a. Hence binary or decimal expansions with ©(|loga|) or more 
digits will be used for all numbers. We will see in the analysis of the algorithm 
below that a precision of 6(|loga|) digits is sufficient. By the Proposition above, 
|loga| < e(logiV). 

Algorithm 

Step 1 By using {2N + 3)1 independent random processes in parallel, generate 
{2N + 3)1 independent samples from the standard Gaussian distribution. Label the 
samples by £(w3,i) and £(a;i,„,i), C{uj2,n,i), ior n = 0, . . . , N, i = 1, . . . ,£. 

Step 2 By using {N + 1)£ independent deterministic parallel processes, compute 
the wavelet coefficients and {ft^KHn) for i = !,...,£ and n = 0,...,N. 

Compute each {f^P ,Hn) = 2^/^ /^'-Tf f^'\u)du - 2^/^ ^^k+H f^'Hu)du, and 
similarly, each {f^f^ , Hn) , by standard numerical integration techniques at a precision 
of 6(1 log a|). 

Step 3 By using deterministic parallel processes, compute for j = 1,. . . ,£ : 

N 

and 

N 

W2{ti,H,N){0J2,i) = Ch J2<fu^^n) ■ C{uj2,n,i) 

n=0 

Here, ujij = {uJi,nA : n = 0, . . . , N), and lo2a = (w2,n,i : n = 0, . . . , N). 

Step 4 By using £ deterministic parallel processes, compute for i = 1, . . . ,i : 

Step 5 By using £ deterministic parallel processes, compute for i = 1, . . . ,^ : 

W{ti , H, N) (wi) = Wi {U , H, N) + W2 {U , H, N) (w2,i) + h {U , H) (wg,,) 



where Wj = (wi,i,W2,j, W3,i)- 
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In step 1, each sample can be generated in constant time, so the parallel algorithm 
also takes constant time. In fact, samples from a standard Gaussian distribution can be 
pre-computed and stored for later use; in that case, step 1 does not need to be counted 
at all. In step 2, each wavelet coefficient can be computed in parallel time 0{\ loga|). 
All these coefficients are computed independently of each other, so the total parallel 
time is 0(| loga|) < 0(log A^). Step 3 carries out 2£ independent additions in parallel 
to find Wj{t„H,N), j = 1,2 and i = For each Wj{t„H,N), N + 1 terms of 

the form {ft^\'Hn )^{'-^j,n,i) are obtained first by independent multiplications, so they 
are obtained (with 0{ \ loga|) digits of precision) in parallel time 0{ \ loga|). The + 1 
terms are then added, which can be done by a parallel algorithm in time 0{\ogN) ■ a, 
where a is the time used to add two items (see e.g., [H]). Addition of two &-bit 
numbers can be done in parallel time 0(log6) using 0{b) processors (see e.g., [11]). 
Thus the parallel time complexity of step 3 is 0(logA^ • log log A^) using O(A^logA^) 
processors. Steps 4 and 5 each use £ independent parallel processes, and each process 
takes time 0{\ loga|) for a precision of | loga| digits. So the time complexity of these 
steps is 0(log A^). Therefore, the total time complexity of the above parallel algorithm 
is 0(logA^ • log log A), using 0{{£ + \ogN) N) processors. 

By the above Proposition we conclude that for a fixed H the parallel time complexity 
is 0(|log(pa)| • log|log(pa)|). 

References 

[1] Arcones, M.A. (1995). On the law of the iterated logarithm for Gaussian processes. Journal of 
Theoretical Probability 8 877-904. 

[2] Bardina, X., JOLIS, M. AND TuDOR, C.A. (2003). Weak convergence to the fractional Brownian 
sheet and other two-parameter Gaussian processes. Statist. Probab. Lett. 65 317-329. 

[3] BlAGlNl, F., Hu, Y., 0KSENDAL, B. AND Zhang, T. (2008). Stocliastic Calculus for Fractional 
Brownian Motion and Applications, Springer. 

[4] Daubechies, I. (1992). Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series in 
Applied Mathematics 61 SIAM. 

[5] Delgado, R. and ,Jolis, M. (2000). Weak approximation for a class of Gaussian process. J. 
Appl. Probab. 37 400-407. 



20 



D. Hong, S. Man, J-C. Birget, D. S. Lun 



[6] Davydov, Y. (1970). The invciriance principle for stationary processes. Tear. Verojatn. Primen. 
15 498-509. 

[7] DZHAPARIDZE, K. AND VAN Zanten, H. (2004). A series expansion of fractional Brownian motion. 
Probability Theory and Related Fields 130 39-55. 

[8] DZHAPARIDZE, K. AND VAN Zanten, H. (2005). Optimality of an explicit series expansion of the 
fractional Brownian sheet, Statistics and Probability Letters 71 295-301. 

[9] Enriquez, N. (2004). A simple construction of the fractional Brownian motion. Stochastic 
Process. Appl. 109 203-223. 

[10] Garzon, J., GOROSTIZA, L.G. AND Leon, J. A. (2009). A strong uniform approximation of 
fractional Brownian motion by means of transport processes, Stochastic Processes and their 
Applications 119 3435-3452. 

[11] Kloeden P.E. and Platen, E. (1999). Numerical Solution of Stochastic Differential Equations, 
Springer. 

[12] KuHN T. AND Linde, W. (2002). Optimal series representation of fractional Brownian sheets. 
Bernoulli 8 669-696. 

[13] Kuo, H-H. (2006). Introduction to Stochastic Integration, Springer. 

[14] Leighton, F. T. (1991). Introduction to Parallel Algorithms and Architectures, Morgan 
Kaufmann. 

[15] Li, Y. and Dai, H. (2011). Approximations of fractional Brownian motion. Bernoulli 17 1195- 
1216. 

[16] Mandelbrot B.B. and van Ness, J.W. (1968). Fractional Brownian motions, fractional noises 
and applications. SIAM Review 10 422-437. 

[17] Meyer, Y., Sellan, F. and Taqqu, M. (1999). Wavelets, generalized white noise and fractional 
integration: the synthesis of fractional Brownian motion. The Journal of Fourier Analysis and 
Applications 5 365-394. 

[18] SzABADOS, T. (2001). Strong approximation of fractional Brownian motion by moving averages 
of simple random walks. Stochastic Processes and their Applications 92 31-60. 



[19] Taqqu, M.S. (1975). Weak convergence to fractional Brownian motion and to the Rosenblatt 
process. Z. Wahrsch. Verw. Gebiete 31 287-302. 



